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The three-body force is indispensable in nuclear energy density functionals which leads to a density depen¬ 
dent two-body term in the Hartree-Fock approach. Usually a single factional power of density dependency 
has been adopted. We consider the possibility of an additional higher-order density dependence in extended 
Skyrme forces. As a result, new extended Skyrme parametertizations based on the SLy4 force are obtained 
and the improvements in descriptions of global nuclei have been demonstrated. The higher-order term can also 
substantially affect nuclear properties in the high density region in general ways. 
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I. INTRODUCTION 


form. 


The nuclear energy density functional theory (EDF) is a 
principal theoretical tool for descriptions of bulk properties 
and dynamics of the entire nuclear landscape. The main is¬ 
sues for nuclear EDF theory are concerns about its accuracy 
and uncertainties. After early Skyrme Hartree-Fock calcula¬ 
tions iQl], the Skyrme interactions have been widely used 
and have been demonstrated to be very powerful JH]. In par¬ 
ticular, the latest UNEDF Skyrme forces based on large 
scale fitting procedures have remarkably improved the accu¬ 
racy. On the other hand, the limits of standard Skyrme forces 
have been reached and new ideas are desirable. 

The Skyrme force is a very low-momentum (or soft) in¬ 
teraction with a second-order momentum dependence and a 
three-body interaction that plays an indispensable role in nu¬ 
clear saturation I0]. In ab initio calculations, it is well known 
that the three-body force and many-body forces are generated 
due to suppressed (or renormalization cutoff) degrees of free¬ 
dom of the two-body potential 1^. Following this picture, 
the remarkable three-body force contribution in soft Skyrme 
forces as well as Gogny forces ||^ is then understandable. 
Recently the three-body force has become an attractive issue 
in ab initio calculations iQiini] and Brueckner-Hartree-Fock 
calculations lfl^[T^ . In Hartree-Fock calculations, the three- 
body interaction can be considered as a density dependent two 
body term with , in which the density dependency power 7 
should be 1. However, the power 7 should be smaller than 2/3 
to get reasonable nuclear matter incompressibility IHIIl. In 
fact, a factional power of density dependency, ranging from 
1/6 to 1, has been adopted in various Skyrme parameteriza- 
tions 10 . The parameterized density dependence is supposed 
to simulate complicated many-body correlations in addition to 
the three-body force and is still an open question. 

In the dilute Fermi systems of repulsive cores, the energy 
density functional can in principle be expanded in powers of 
kp (or p'/^). The EDF expansion coefficients were calcu¬ 
lated in a series of papers in 1957, by Huang-Yang ifl^ . Lee- 
Yang ifl^ . Martin and De Dominicis ifl^ . which takes the 
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in which a denotes a finite radius of particles. These high 
order terms reflect the increased kinetic energy by consid¬ 
ering the finite radius of particles. In the Local Density 
Approximation for dilute Fermi gases, its EDF form is the 
same as Skyrme EDF except for different density dependent 
terms lEoll . The EDF of dilute Fermi gases involves terms 
with p^/®, p®/®, etc. While the standard Skyrme EDF involves 
a single density dependent term with p^+^. Actually, the p^/® 
term (corresponding to 7 = 1 / 3 ) can be mimicked by an effec¬ 
tive three-body interaction In contrast to the dilute Fermi 
systems, we note that generally the finite radii (volume ef¬ 
fects) of nucleons have not been explicitly considered, neither 
in realistic nor phenomenological nuclear forces. 

Our objective of this work is to consider the possibility of 
an additional higher order density dependence on top of the 
standard Skyrme force. In deeded, in the standard Skyrme 
force, a single density dependent term might be too simplis¬ 
tic. It is also of interdisciplinary interests to study the role 
of higher order density dependencies and the EDF connection 
between nuclei and dilute Fermi gases. We speculate that the 
high-order density dependencies can reflect short-range corre¬ 
lations beyond standard Skyrme forces, and impact in partic¬ 
ular the high-density behaviors and isovector properties. 

The extensions of Skyrme forces have been developed by 
many efforts. For example, the Brussels forces by considering 
the density dependence in the momentum-dependent terms 
have been very successful The effective pseudopoten¬ 
tial of Skyrme-type EDF has been derived up to the sixth or¬ 
der . On the other hand, the connection between nuclear 
EDF parameters and ab initio calculations has not yet been 
clear, although progresses have been achieved in the Density 
Matrix Expansion and the Effective Field Theory ll2^ . 

The higher-order density dependent terms have been con¬ 
sidered in Refs, 1^ in analogy to those for the dilute 
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Fermi gases, but with limited applications. In this work, we 
are not intended to refit a fully new Skyrme parameter set from 
scratch but to improve the existed successful Skyrme param¬ 
eters with an additional higher-order density dependent term. 
To this end, we concentrate on the momentum-independent 
terms, and manipulate the to, to parameters and the extended 
term with t^E- Indeed, to, to and toE are directly entan¬ 
gled in the s-wave channel although they are also connected 
with other parameters. For example, the correlations between 
parameters have been demonstrated by the Bayesian analy¬ 
sis 12^ . We choose to refit the momentum-independent terms 
based on the SLy4 force and keep other parameters un¬ 
changed. The SLy4 force is an ideal start because it has been 
widely used and has been the references for other parame- 
terizations such as UNEDF forces. Furthermore, SLy4 has 
a density dependence of 7 = 1/6 which is relatively small and 
may leave room for higher-order density dependencies. Other 
successful Skyrme forces such as SkM* with 7 = 1/6 may 
also be improved with higher-order density dependencies. Or 
say we consider an additional higher-order density dependent 
term for complementary to SLy4 in a perturbative manner. It is 
turned out that the density dependent terms we employed are 
different from Refs. j^l^ . We will then evaluate the per¬ 
formance of the extended new parameterizations compared to 
the original SLy4. 


11. THEORETICAL METHODS 


We begin with the standard Skyrme interaction including 
the two-body and three-body terms || 2 l, which takes the form. 


hed as. 
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In this case, there will be 2 more additional parameters toE 
and xoE- In the SLy4 force the power factor 7 takes 
1/6 and then we consider the next higher order power of 
1/2 = 1/6-1- 1/3. Actually such an extension is straightfor¬ 
ward and can be easily implemented. We noticed that similar 
extensions exactly according to the density dependent terms 
for dilute Fermi gases have been discussed in the case of nu¬ 
clear matter properties 12^ . Later it has been investigated in 
finite nuclei but negative coefficients are obtained for high- 
order terms . The density dependent terms we employed 
are different from ll^ | 2 ^ since our studies are based on 
SLy4, although we have similar motivations. It is useful to 
explore different combinations of density dependencies. Con¬ 
sidering recent UNEDF forces with 7 around 1/3 US] that 
is between 1/6 and 1 / 2 , the two density dependent terms we 
employed should be reasonable. 

Next we refit the extended Skyrme parameters for finite nu¬ 
clei with the Simulated Annealing Method j^ . In this work, 
we only refit the momentum-independent parameters, to, to, 
toE, xq, Xo, Xoe, and keep others ti, t 2 , xi, x^. Wo un¬ 
changed based on the SLy4 force. In this way the influences 
of the extended higher-order density dependent term can be 
clearly illustrated, to avoid the influences due to terms which 
involve momentum dependencies. 

Our fitting procedure is similar to SLy4 as described in 
Ref. ||25] . Briefly, we minimize the quantity. 


i<j i<j<k 
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The three-body interaction in the Hartree-Fock calculations 
can be transformed into a density dependent two-body inter¬ 
action with . Correspondingly the power 7 should be 1 but 
usually it is adjusted to be less than 1 in modern Skyrme forces 
to get reasonable incompressibilities. 

= ^ 3(1 + xoPa)5{ri - rj)S{rj - r^) 

( 2 )' 1 

^ Wj = gi3(l + XoPcT)p(Rr<^(ri-rj) 

In Eq.© and Eq.® , t^, Xi and Wo are the parameters of the 
Skyrme interaction, and R = (r^ -b 'rj)l2. 

In our extension, by including an additional higher-order 
density dependent term, the density dependent term is modi- 
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where Coo is the average energy per nucleon at the satura¬ 
tion point; and Euvi 4 +uvii{i) energies of the neu¬ 
tron matter from Skyrme forces and calculations of Wiringa 
et al. lIJIl] at different density points; Bt and Rc denote total 
binding energies and charge radii of selected nuclei respec¬ 
tively. In Eq.®, the htting uncertainties AEi are not a con¬ 
stant but increase as densities increase M- The nuclear mat¬ 
ter properties at the saturation point are important inputs for 
constraining parameters. In the extended Skyrme functional, 
we have to slightly modify the equations of the pressure P, 
the incompressibility coefficient K, the symmetry energy co¬ 
efficient as- The saturation density of 0.1595 fm“^, the sym¬ 
metry energy coefficient as of 32 MeV, and the effective mass 
m*/m of 0.7 have not been adjusted 11^ . Therefore, we ad¬ 
just 6 parameters and have 4 free parameters Coo, to, xo and 
Xoe to fit after taking into account the equations of nuclear 
matter properties. Then Skyrme parameters to, xo, toE can be 
determined through these expressions. 

Our calculations are based on the axially-symmetric 
Hartree-Eock-iBCS approach and the Simulated Annealing 
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FIG. 1. (Color online) The optimization of Skyrme parameter sets 
by adopting different fitting groups with the Simulated Annealing 
Method. The obtained values are displayed as a function of 1/T, 
where T denotes the temperature in SAM. 



Density (fnn‘^) 

FIG. 2. (Color online) Energies per particle in nuclear matter ob¬ 
tained by different Skyrme parameter sets (see TableUll as a function 
of density. Yp=0 denotes for the pure neutron matter and 1^=0.5 
denotes for the symmetric nuclear matter. 


Method (SAM). The Hartree-Fockn-BCS equation is solved 
within axially-symmetric coordinate spaces using the Skyax 
solver 1^ . In this case deformed nuclei can be included in 
the fitting procedure compared to the original SLy4 force. The 
calculations are performed in a 2D box of 30 fm and the uni¬ 
form lattice space is taken as 0.5 fm. The mixed pairing inter¬ 
action 13^ has been adopted and the pairing strength is taken 
as Vp^n=500 MeV fm“^ for protons and neutrons. SAM 
is a general purpose algorithm for multi-parameter fittings, 
based on the Monte Carlo iterative solution method. SAM 
starts with a high temperature and randomly searches the min¬ 
imum, which accepts a worse solution than the current one 
with a certain probability. As the temperature cools down, the 
searching space becomes smaller and the optimization is then 
realized. The SAM has already been adopted in the fitting of 
full Skyrme parameters . 


III. RESULTS AND DISCUSSIONS 



FIG. 3. (Color online) Incompressibility coefficients K as a function 
of nuclear matter density, corresponding to different Skyrme param¬ 
eter sets in TableU 


We refit the extended Skyrme force based on SLy4 with 
different groups of selected nuclei. In particular, we select 4 
groups aiming at different mass regions, 

Groupl: ^°Ca,‘*®Ca,^®Ni, ^°®Pb, for magic nuclei 

as adopted in the SLy4 fitting; 

Groupl; i6O,36Mg,40Ca,4®Ca,^°Cr,56Ni,7®Ni,i00Zr, 

for light nuclei with 8 ^ ^ 50; 

Group3: 120Sn,132Sn,l™Gd,l™Hf,198pb,208pb,236u^252] 
^™Sg, for heavy nuclei with Z ^ 50; 

Group4: i6O,36Mg,40Ca,4®Ca,^°Cr,56Ni,7®Ni,i00Zr, 

i20Sn,i32Sn,i60Gd,i™Hf,i98pb,208pb,236u,252Fm, 260Sg, 

for the global mass region. 

The exp erimental binding energies of these nuclei are taken 
from On]. The charge radii of 5 magic nuclei are taken 
from 11 ^ . 


Figure [T] illustrates the obtained of different SAM fit¬ 
tings as the temperature cools down. In the fitting procedure, 
the initial temperature is 2 and decreases by a factor of 0.5 
each step. The initial parameters are given as random values in 
setted ranges. The initial ranges of X 3 and x^e are from —2 to 
2. The range of Coo is 16±0.15 MeV. The range of ts is from 
8000 to 13777. At the beginning, the is very large and be¬ 
comes much smaller as the temperature lower than 0.125. One 
. full fitting procedure takes about 200 iterations for the fitting 
of light nuclei and about 350 iterations for the global fitting, 
which takes 2'^3 days. After that, we polish these parameters 
by a second fitting with obtained parameters as initial values. 

TableU] lists the obtained different parameterizations of ex¬ 
tended SLy4 forces compared to the original SLy4 force. Gen¬ 
erally, the obtained x^ values of charge radii and binding en¬ 
ergies are reduced significantly, compared to the values of 
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TABLE I. Refits of SLy4 parameters with extend density dependent terms, compared to the original SLy4 force. The different parameter sets 
are obtained by fitting different groups of nuclei, see text for details. We also refit SLy4 without the extended term as labeled by SLy4'. Other 
parameters in SLy4 not listed here are not adjusted. The corresponding minimized -x^ of binding energies (Xb) and charge radii (Xb) of SLy4 
are given in the brackets. 


Parameters 

SLy4 

SLy4'(global) 

Group 1 (magic) 

Group2(light) 

Group3(heavy) 

Group4(global) 

fo(MeV • fm^) 

-2488.91 

-2493.536 

-2132.57 

-2106.69 

-2310.79 

-2319.15 

f3(MeV-fm®(^+H^) 

13777.0 

13809.324 

9366.12 

9036.30 

11549.45 

11660.70 

f3B(MeV • 

0 

0 

2756.97 

2970.72 

1410.64 

1334.88 

Xo 

0.834 

0.931 

0.968 

0.988 

0.972 

0.863 

X3 

1.354 

1.496 

1.922 

1.975 

1.715 

1.509 

X3E 

0 

0 

0.375 

0.530 

0.324 

0.622 

Coo (MeV) 

-15.972 

-15.963 

-16.063 

-16.085 

-16.055 

-16.039 

Ko, (MeV) 

229.9 

230.31 

247.33 

248.84 

239.19 

238.58 

Xr 


10.35 

2.18(4.04) 

5.00(10.04) 

0.39(0.62) 

7.4(10.66) 

Xb 


31.67 

0.83(2.46) 

3.83(6.82) 

19.26(78.28) 

25.92(84.87) 

Xsum 


42.84 

3.67(6.91) 

9.79(17.306) 

20.52(79.21) 

33.82(95.94) 


TABLE 11. Energies (in MeV) of neutron droplets in a tkj=10 MeV 
spherical trap calculated with SLy4 and the extended SLy4 (light) 
parameters, compared to ab-initio results The extended SLy4 
parameters are given in Table |I] 


N 

SLy4 

Extend (light) 

ab initio 

8 

120.81 

124.03 

135.4(2)^ 

10 

169.54 

173.47 

187.1(6)“ 

12 

217.04 

221.56 

237.1(8)“ 

14 

262.69 

267.59 

286(1)“ 

16 

309.28 

314.28 

334(1)“ 

18 

364.83 

370.87 

386(3)“ 

20 

418.56 

425.50 

432.8(5)^ 

28 

656.18 

663.84 

681(1)'’ 

40 

1061.23 

1070.97 

1058(2)*’ 

50 

1420.19 

1428.28 

1449(3)'’ 


^ No-core shell model calculations 
Couple-cluster calculations 


SLy4 shown in brackets. The incompressibility values 
and energies Coo at the saturation density are also listed. We 
can see that the existed parameters to, are reduced. The 
obtained t^E parameter of the higher-order term ranges from 
1334.88 to 2970.72. We noticed that in Refs, ll^ [^ . with 
decreasing momentum cutoffs in second-order Hartree-Fock 
calculations, the adjusted parameters to, to decrease and the 
density dependency power increases. Our results with higher 
order terms are similar to that trend. In the fittings of heavy 
and global nuclei, the parameters of the higher-order term are 
relatively small. While in the fitting of light nuclei, the param¬ 
eter of the higher-order term is the largest, which is mainly re¬ 
sponsible for improving the descriptions of charge radii. This 
indicates that complex density dependencies are required for 
describing light nuclei. In our tests, by reducing the fitting 
weights on charge radii, the obtained higher-order term de¬ 
creases. In the heavy mass region, the obtained new extended 
parameterizations still produce large deviations in binding en¬ 
ergies. We also refit the SLy4 with the global group, as labeled 
by SLy4'. The resulted is significantly reduced except for 



FIG. 4. (Color online) Symmetry energies as a function of nuclear 
matter density, corresponding to different Skyrme parameter sets in 
TableH 


the descriptions of charge radii. This demonstrated that the 
extended higher-order density dependent term is important for 
descriptions of nuclear charge radii. 

Figure|2]displays the performances of extended parameter¬ 
izations for energies of the neutron matter and the symmetric 
nuclear matter. For the neutron matter, the energies given by 
UV14H-UVII are also shown. Since our fittings have in¬ 
cluded the UV14H-UVII results, we see energies are close to 
the original SLy4. In the case of neutron matter, energies of 
all the extended parameterizations are slightly decreased in 
the high density region compared to the original SLy4. While 
in the case of symmetric nuclear matter, energies are slightly 
increased. Note that large uncertainties in the high density re¬ 
gion have been taken into account in the fitting procedures. 
The energies of neutron droplets in a fiuj=lQ MeV trap have 
also been calculated, as displayed in Table In] The ab initio 
calculations were done by the no-core shell model and the 
couple-cluster model with the Chiral two-body (N3LO) plus 





















5 



0 20 40 60 80 100 120 140 160 

Neutron number 


FIG. 5. (Color online) Binding energy differences between theory 
and experiments (H for 603 even-even nuclei. The results are ob¬ 
tained by Hartree-Fock-l-BCS calculations with (a) the SLy4, (b) the 
refitted SLy4' and (c) the extended SLy4 (global), respectively. See 
TableUfor the parameter sets. 


three-body (N2LO) forces 1^ . Generally, the SLy4 force un¬ 
derestimates the energies compared to the Chiral force. While 
the extended SLy4 force can bring more binding energies. 
This is consistent with the situation of neutron matter ener¬ 
gies as shown in Figl2l 

To study the behaviors of extended SLy4 forces for nuclear 
matter properties, the incompressibility and the symmetry en¬ 
ergy are displayed in Fig[3]and FiglUrespectively. In FiglJl the 
extended forces all generate larger incompressibility values at 
high densities compared to the SLy4 force. The parameters 
for light nuclei with the largest higher order term produce the 
highest incompressibility. This causes very little influences in 
the low density region. It is understandable that the increased 
incompressibility is due to the repulsive higher-order density 
dependence. In FigH] the symmetry energies as a function of 
densities are displayed. It can be seen that the symmetry en¬ 
ergies decrease in the high density region by considering the 
higher order density dependent term. Again the parameters 
for light nuclei produce the softest symmetry energy. There 
has been some favorable arguments about the soft symmetry 
energy 1 ^ . in contrast to the observation of a neutron star 
with two solar mass. While heavy ion collision experiments 
can not give clear constraints about the symmetry energy in 
dense nuclear matter S. In fact, by studying the expres¬ 
sions of nuclear matter properties, we note that the behav¬ 
iors of the extended higher-order density dependence shown 


in Figs. (I34 |i are general. We conclude that the inclusion of a 
higher-order density dependent term would impact the equa¬ 
tion of state and isovector properties in the high density region 
(supra-saturation density), which are important for studying 
neutron star structures. 

Finally, the global calculations of binding energies of 603 
even-even nuclei are shown in FiglS] The calculations are 
done with the axially-symmetric Hartree-Fockn-BCS solver 
Skyax, and with the original SLy4, the refitted SLy4' and 
the extended SLy4 (global) parameter sets. As we can see that 
the results of the extended parameterization are similar to the 
SLy4' results. The original SLy4 remarkably underestimates 
the binding energies of heavy and superheavy nuclei. The ex¬ 
tended SLy4 and SLy4' have improved the global descriptions 
by including some heavy nuclei in the fitting procedure. For 
all the nuclei, the resulted root-mean-square (rms) deviations 
of the extended SLy4, SLy4' and SLy4 are 2.30 MeV, 2.94 
MeV and 4.37 MeV, respectively. In the region of Z < 82, 
the rms deviations of the extended SLy4, SLy4' and SLy4 are 
2.00 MeV, 2.71 MeV, 3.18 MeV, respectively. In the region of 
Z ^ 50, the rms deviations of the extended SLy4, SLy4' and 
SLy4 are 2.15 MeV, 2.88 MeV, 1.97 MeV, respectively. There¬ 
fore we can say that the extended SLy4 force has globally im¬ 
proved the descriptions of binding energies, compared to the 
refitted SLy4' and the original SLy4. Note that the extended 
global parameterization set has a relatively small higher-order 
density dependent term, as displayed in Table U We see that 
SLy4, SLy4' and the extended SLy4 are not able to describe 
the binding energies of magic nuclei and deformed nuclei si¬ 
multaneously. Or say the shell effects are overestimated as 
they are based on a small effective mass of 0.7 ll^ . Presently 
the effective mass and the related ti, (2 parameters have not 
been adjusted in the extended SLy4 forces. Indeed, UNEDF 
Skyrme forces with large effective masses around 1 are very 
successful in describing binding energies of the whole nuclear 
landscape 10]. On the other hand, the Brussels Skyrme force 
with an effective mass of 0.8 and additional collective cor¬ 
rections for deformations also works very successful lli^ . In 
addition, SLy4 is unable to describe the surface tension com¬ 
pared to calculations with microscopic center-of-mass correc¬ 
tions lEH. Actually, our main goal in this work is to study the 
influences of the extended higher-order density dependence, 
by refitting only the momentum-independent parameters. The 
global performance of the extended SLy4 is very encouraging. 
To further improve the overall performance, the adjustments 
of all the Skyrme parameters including the extended term by 
choosing properly physics inputs will be our next step. 


IV. SUMMARY 

In summary, we studied the extension of Skyrme forces by 
including a higher-order density dependent term, according to 
our speculation that a single density dependent term might be 
too simplistic. This is in analogy to high order density de¬ 
pendent terms in the EDF for dilute Eermi gases. We studied 
the influences by adjusting only the momentum-independent 
parameters based on the SLy4 force, with two density depen- 
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dent terms of and The obtained new extended pa- 
rameterizations indicate that the strengthes of the higher or¬ 
der term p^/^ are dependent on different fitting regions. We 
have seen that the higher-order term can improve descriptions 
of binding energies of global nuclei and neutron droplets as 
expected. We also noticed that the higher-order term is im¬ 
portant for descriptions of charge radii. The extended param¬ 
eterization set obtained by fitting light nuclei has the largest 
higher order term, indicating complex density dependencies 
are required for light nuclei. Furthermore, we demonstrated 
that the extended density dependence can impact nuclear mat¬ 
ter properties particularly in the supra-saturation density re¬ 
gion with general behaviors, although the extrapolations of 
Skyrme forces to the high density region still have large un¬ 
certainties. In conclusion, our studies have provided some in¬ 
sights and opportunities for future developments of Skyrme 


energy density functionals. 
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